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The dynamics of a set of rods bouncing on a vertically vibrated plate is investigated using ex- 
periments, simulations, and theoretical analysis. The experiments and simulations are performed 
within an annulus to impose periodic boundary conditions. Rods tilted with respect to the vertical 
are observed to spontaneously develop a horizontal velocity depending on the acceleration of the 
plate. For high plate acceleration, the rods are observed to always move in the direction of tilt. 
However, the rods are also observed to move opposite to direction of tilt for a small range of plate 
acceleration and rod tilt. A phase diagram of the observed motion is presented as a function of 
plate acceleration and the tilt of the rods which is varied by changing the number of rods inside 
the annulus. Next we introduce a novel molecular dynamics method to simulate the dynamics of 
the rods using the dimensions and dissipation parameters from the experiments. We reproduce the 
observed horizontal rod speeds as a function of rod tilt and plate acceleration in the simulations. 
By decreasing the friction between the rods and the base plate to zero in the simulation, we identify 
the friction during the collision as the crucial ingredient for occurrence of the horizontal motion. 
Guided by the data from the experiments and the simulations, we construct a mechanical model for 
the dynamics of the rods in the limit of thin rods. The starting point of the analysis is the collision 
of a single rod with an oscillating plate. Three friction regimes are identified: slide, slip-stick, and 
slip reversal. A formula is derived for the observed horizontal velocity as a function of tilt angle. 
Good agreement for the horizontal velocity as a function of rod tilt and plate acceleration is found 
between experiments, simulations and theory. 



I. INTRODUCTION 

Granular materials come in all shapes and sizes. Ide- 
alized spherical particles, have been typically used to un- 
ravel the fascinating properties displayed by granular ma- 
terials. In vibrated granular systems, periodic pattern 
formation, cluster formation, and complex size separa- 
tion have been observed. However, anisotropic grains are 
nearly as numerous, and experience with thermal systems 
teaches us that shape matters. Only a handful of investi- 
gations have studied the impact of anisotropy on granular 
systems, but the ones that have been accomplished point 
to a rich phenomenology. 

For example, theoretical and numerical study of a low- 
density gas of hard inelastic needles Q shows two dis- 
tinct regimes of cooling related to the different scaling of 
rotational and translational energy dissipation. In com- 
paction experiments with granular rods |3] granular rods 
vibrated in a tall narrow container were observed to form 
ordered stacks i.e. a smectic phase similar to that found 
in thermal systems, with the additional novelty that the 
rods align with the gravitation field. More recently, self- 
organization of vortices was observed to occur when a 
shallow bed of granular rods was vibrated Q ■ It was fur- 
ther shown, that the tilt of the rod and vertical vibration 
was important to the occurrence of the novel dynam- 
ics. While a phenomenological model of formation and 
growth of the vortices has been proposed 0| , a detailed 
understanding of why the rods move horizontally on a 
vertically vibrated plate was not reached. 

The collective motion of vibrated anisotropic grains 
is of considerable interest as an example of spontaneous 
ratchet formation in a non-equilibrium dissipative sys- 



tem. Transport of thermal particles in systems with mi- 
croscopically asymmetric potential has been studied in 
a number of recent publications JJ. In Ref. 0, ratchet 
transport was demonstrated for spherical grains on a ver- 
tically vibrated asymmetric saw-tooth profile. In Ref. 0, 
the transport of elongated grains on a vertically vibrated 
ratchet-shaped plate has been studied numerically. How- 
ever, as follows from results of 0, 0, and further de- 
scribed in this paper, the transport of anisotropic grains 
may occur even without externally imposed microscopic 
asymmetry of forcing, in which case the direction of mo- 
tion is chosen as a result of spontaneous symmetry break- 
ing. 

In this paper, we apply experimental and numerical 
tools to a system of rods in a vibrated annular container 
to elucidate the development of coherent horizontal dy- 
namics in anisotropic systems. This geometry was specif- 
ically chosen to simplify the phenomenology in order to 
focus on the mechanisms for the observed dynamics. Our 
theoretical model is developed for even simpler quasi-two- 
dimensional geometry where all rods are confined to a 
vertical plane, and periodic boundary conditions are im- 
posed. The theory is based on a detailed description of 
frictional collisions between rods and the vibrating plate 
which makes use of the assumption of constant kinematic 
restitution coefficient. While the issues of restitution co- 
efficient in application to frictional impact of asymmetric 
bodies are still debated in the literature (see, for example, 
'^), we show that even this simplest model agrees very 
well with soft particle molecular dynamics simulations of 
individual collisions. To describe the collective motion of 
the rods, we take advantage of the observation that in the 
regime of stationary translation the mean horizontal mo- 



mentum transfer due to the collision with bottom plate 
is zero, and assume that this condition holds for a typi- 
cal collision. Our simulations show that during the flight 
the angular momentum of a single rod is transfered to 
other rods, so the angular velocity of the rod at the end 
of the flight becomes small and can be neglected. Fur- 
thermore, based on our numerical simulations we make 
the assumption that the vertical velocity of the center of 
mass (CM) just before the collision is equal to the CM 
velocity in the beginning of the flight. These assump- 
tions allow us to find the mean translation CM velocity 
in a closed form. We show that this theory captures the 
essential mechanisms of the transport of tilted rods on a 
vibrating plate. 

The paper is organized as follows. First, we introduce 
the experimental system and report the observed dynam- 
ics as a function of the control parameters such as plate 
acceleration and tilt of rods. Next we discuss the molecu- 
lar dynamics simulations corresponding to experimental 
parameters and compare the results with the experimen- 
tal data. Then using the data as a guide, we construct 
a theoretical model for the occurrence of spontaneous 
horizontal motion, and its dependence on control param- 
eters, and compare the results with the data. Finally, we 
discuss the general implications of this study. 



II. EXPERIMENTS 

The experimental system consists of an annular con- 
tainer with an inner diameter of 7.28 cm and outer di- 
ameter of 9.45 cm, and is similar to that used in Ref. 
An image is shown in Fig. . The sides are composed 
of clear acrylic and the base plate is made of Aluminum. 
The container is attached to an electromagnetic shaker 
through a linear bearing which allows only vertical mo- 
tion. A frequency generator along with a lock-in amplifier 
is used to excite the system with a fixed frequency and 
amplitude. The data reported here for frequency / = 60 
Hz, and we note that qualitatively similar behavior is ob- 
served when the frequency is varied between 50 and 100 
Hz. The acceleration of the container is measured with 
an accelerometer and is reported in terms of F, the mea- 
sured peak acceleration divided by the acceleration due 
to gravity. The tilt to the rods is characterized by the 
angle (f) with respect to the vertical axis. 

The rods used for the experiments are cylindrical with 
a diameter of 0.635 cm and length 5.08 cm. One of the 
ends is semi-spherical with the radius equal to the radius 
of the rod, and the other end is flat. The rods are made 
of a Delrin and Teflon composite, and the measured dissi- 
pation coefficients are as follows. The coefficient of static 
and kinetic friction between the rod and the base plate 
is 0.36 and 0.25 respectively. The coefficient of restitu- 
tion is obtained by measuring the kinetic energy of the 
rod just before and after a collision with a stationary 
base plate. For nearly normal incidence (0 < 5'^), the 
coefficient of restitution is approximately 0.8 ±0.1 (see 
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FIG. 1: (a) Annular geometry used in the experiments, (b) 
The horizontal and vertical position of the rod end at the 
bottom as a function of time. The position of the oscillating 
bottom plate is also plotted for reference. The three plots are 
shifted for clarity. (F = 3.3, = 50, = 13") 



also ))\ .) When a single rod is placed inside the annulus, 
the maximum tilt angle that it can have is 53*^ due the 
curvature of the annulus. By increasing the total number 
of rods N in the annulus from = 23 to 56, the tilt (f> 
can be varied from 53° to O''. 

FigureQ^b) shows the typical motion of the bottom tip 
of a rod as a function of time. The data is obtained by 
using a high-speed Kodak digital camera with a frame 
rate of 1000 per second, and tracking the end of the rod 
with appropriate use of lighting through the transparent 
side walls. The vertical position of the tip z is observed 
to oscillate, as the rod bounces on the vibrating plate 
[also shown in the Fig. ^b)]. The flight time is observed 
to vary and although the rod appears to almost always 
hit the plate on its upstroke, a distribution of phases is 
observed. We will discuss this issue further after intro- 
ducing the molecular dynamics simulations in a later sec- 
tion. On the other hand, the horizontal position of the tip 
is observed to increase approximately linearly although 
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some oscillatory motion is also seen be superposed. The 
slope of the x-position gives the average horizontal veloc- 
ity which is consistent with dividing the average circum- 
ference of the annulus by the amount of time taken by 
the rods to go around once. This second method is used 
to report observed average horizontal velocities. 

The measured horizontal speed Cx as a function of the 
tilt of the rods cj) for a fixed F = 3.3 is shown in Fig.|2Ia). 
Cx is observed to increase from zero to a peak value and 
then decrease. The rods are observed to always move 
in the same horizontal direction as the tilt. It is to be 
noted that the error in determining Cx arises from run to 
run variability due to slight differences in packing inside 
the annulus rather than from the actual measurement of 
the velocity itself. In this case, seven separate realiza- 
tions were used to arrive at the average value of Cx for a 
particular number of rods. 

By varying F for a few values of 0, we obtain its impact 
on observed horizontal velocities [see Fig. Efb) .] Positive 
horizontal velocity is observed to commence only above a 
finite F ~ 1.6 is reached. Below F 1.6, no net horizon- 
tal velocity occurs except for a narrow parameter range 
where a horizontal motion in a direction opposite to the 
tilt is observed. As shown in the inset to Fig.|2Ic), this 
reverse horizontal motion is observed to be two orders of 
magnitude slower than the forward motion. 

A phase diagram of the various kinds of observed mo- 
tion is shown in Fig. I3d) as a function of F and cj). For- 
ward motion indicates parameters for which horizontal 
motion occurs in the direction of tilt, and reverse mo- 
tion, indicates when the motion is in the opposite direc- 
tion. The reverse motion is observed to occur only for 
a narrow range of parameters. The horizontal motion 
is predominately along the direction of tilt provided a 
minimum acceleration for the container is achieved. 

To study the impact of the shape of the rod tip, we 
also performed experiments by flipping the rods so that 
the flat end is at the bottom. Figure [2Ic) compares the 
measured Cx as a function of F for rods with rounded and 
flat ends but otherwise under identical conditions. When 
the flat end interacts with the bottom plate, the observed 
horizontal velocities are lower, the minimum rod tilt re- 
quired to obtain horizontal motion is 8°, but otherwise 
the qualitative phenomena remains the same. 
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III. NUMERICAL SIMULATIONS 

We performed a series of numerical simulations of the 
bouncing rods on a vibrated plate. The rods are mod- 
eled as spherocylinders of diameter d, length mass m, 
and moment of inertia /. A rod has three translational 
and two rotational degrees of freedom, the rotation of 
a rod around its own axis is neglected. Our numerical 
algorithm is based on the "soft spheres" molecular dy- 
namics technique TO]. The interaction forces between 
colliding spherocylinders are calculated via the interac- 
tion between viscoelastic virtual spheres of diameter d 



FIG. 2: (a) The average horizontal velocity Cx of the rods as 
a function of tilt angle 0. (b) Cx versus F for three different 
tilt angles, (c) The effect of rod end shape on measured c^. 
A flat end rod is observed to move slowly compared to a rod 
with a round end. Horizontal motion in a direction opposite 
to the tilt is observed over a small range of low F. {cj> — 28° ± 
2°) Inset: Cx measured between F = 1 and 2 is replotted to 
clarify the small magnitude of the reverse motion, (d) A phase 
diagram denoting the kinds of horizontal velocity observed 
relative to the direction of the tilt. The rods are observed to 
move in the direction of the tilt under most conditions. 
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FIG. 3: Two types of geometry used in simulations: (a) - an- 
nulus geometry, (b) - quasi-2D geometry with periodic bound- 
ary conditions along x where the motion of rods is constrained 
to the X — z plane. 



centered at the closest points between the axes of sphe- 
rocylinders, so that the cylinders are in contact whenever 
virtual spheres are. The normal forces between virtual 
spheres are computed using Hertzian model and the tan- 
gential frictional forces are computed by the Cundall- 
Struck algorithm. They lead to the kinematic restitution 
coefficient of about 0.65-0.7 slightly varying with impact 
angle. In most cases we used equal friction coefficients 
A'rr = Mrb = Mrui = 0.3 for all sliding surfaces (rod- 
rod, rod-bottom and rod-side wall) and ignored the dif- 
ference between dynamic and static friction coefficients 
while comparing to experiments. The forces arising from 
the interaction of virtual spheres are then applied to the 
rods (see Appendix A for details). The motion of rods 
was obtained by integrating the Newton's equations with 
the forces and torques produced by interactions of a rod 
with all the neighboring rods, walls of the container, and 
by gravity. 

Figure|2|illustrates two configurations employed in nu- 
merical simulations. An annular geometry [Fig.|3fa)] was 
used to match closely the experimental setup. However, 
to separate the effects of side wall friction and oblique 
collisions among the rods due to annulus curvature, we 
also studied the quasi-2D geometry of Fig. Olb) in which 
axes of all rods are constrained to the x — z plane. This 
geometry is more amenable to the theoretical analysis 
and was used for comparison with our theoretical predic- 
tions. 

In simulations we observed robust drift of the rods 
in the direction of inclination in agreement with exper- 
iments ■ Figure 01a) shows the translation velocity 
as a function of the inclination angle for the annulus ge- 
ometry and parameters / = 60 Hz and F = 3.3 used 
in experiments. The translation velocity grows linearly 
for small 0, reaches maximum at </) « 18° and also for 
(j) « 35° after which it decays to zero at large 0. For in- 
termediate inclinations, 18° < cj) < 35°, in most cases we 
observed noticeable slowdown which however is not typi- 
cally observed in experiment. We will discuss the source 
of this discrepancy below. 

Figure^fb) shows the dependence of the average trans- 
lational velocity of rods on the acceleration of the con- 
tainer at fixed frequency for a number of inclinations. As 
in experiments the motion starts above a (slightly lower) 
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FIG. 4: The results of simulations in the annular geometry: 
(a) - average translational velocity of rods as a function of 
tilt at / = 60Hz,T — 3.3; (b) - idem as a function of F for a 
number of tilt angles; (c) - idem and average tilt as a function 
of the coefficient of friction between the rods and the walls; 
(d) - idem and average tilt as a function of the coefficient of 
friction between the rods and the bottom. 
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FIG. 5: The results of simulations in the annular geome- 
try illustrating formation of the non-uniformity in the tilt for 
/ = 60Hz, r — 3.3. Shown are (normalized) horizontal pro- 
jections of the rods' directors for two numbers of rods (a) 
N = 55, (b) iV = 44, the mean tihs are (j> ^ 9° and <p k: 22° 
respectively. 



threshold F > 1.5. Above F « 2.0 this velocity grows 
roughly linearly with the acceleration. We explain the 
presence of the threshold by the friction with side walls. 
Indeed, in the experiments in quasi-2D geometry there is 
no threshold and c^; oc F — 1 right down to F = 1.0. Fig- 
ure ^fc) depicts explicitly how the average translational 
velocity and average tilt depend on the coefficient of fric- 
tion with the sidewall /i^tu, for a fixed number of rods. As 
one could expect this dependence shows monotonic decay 
when fj,rw is increased; less expected is more than ten- 
fold difference in velocities for large and small firw which 
probably explains systematically slower translational ve- 
locities in simulations and underlines an importance of 
the proper account of friction with walls. 

Next, to eliminate the effects of the curved side walls 
and out-of-plane rod-rod collisions, we simulated the col- 
lective rod motion in a quasi-2D geometry with periodic 
boundary conditions along the direction of the rod tilt 
[Figure |3^b)]. All rods are confined strictly to the x — z 
plane, and the interaction with side walls was ignored, 
while friction with the vibrating bottom plate and among 
the rods was taken into account. We used a fixed number 
of rods {N = 40) and varied the length of the container, 
thereby changing the tilt angle (j) [see Figure OIc)] . The 
relation between the length and the tilt is well described 
by a simple formula cos0 — dN/L. 

Figure EI a-c) the mean values of rod velocities before 
and after impacts. As seen from Fig. |H{a), the angular 
velocity before the collision is rather small (apparently, 
it decays after inelastic collisions with other rods during 
flight). Figures ini[b,c) shows the horizontal and vertical 
CM velocities of rods just before and after collision in the 
laboratory frame as a function of (f>. As seen from these 
plots, the pre- and post-coUisional velocities are close to 
each other. The mean velocity of the plate at the moment 
of collision is shown in Fig. EJd). The velocity is only 
weakly dependent on (p and is close to Vo/2. In Fig.|B|[a) 
one can also see a noticeable variation of the horizontal 
translation velocity near (j) = 30° which however is not 
accompanied by a sharp drop at « 20° observed in 3D 
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FIG. 6: The mean CM velocities of rods before (black) and 
after (red) collision as a function of inclination angle </!> in 
quasi-2D rmmerical simulations. The parameters are I = 9.5d, 
r = 3.3, fi = 0.3. For comparison with experiments the 
velocities are plotted in dimensional units for experimental 
values of d = 0.6cm and / — 60Hz: (a) angular velocities uj 
and uj'; (b) vertical CM velocities — Cz and c'^; (c) horizontal 
CM velocities and c^; (d) the mean plate velocity at the 
time of collision. 



We conclude that there are two different mechanisms 
which independently contribute to the slowdown of rods 
at intermediate values of tilt angle. The first mecha- 
nism that only operates in 3D geometry is related to the 
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FIG. 7: Typical trajectories of the tip of a rod in numerical 
simulations for two different tilt angles (\> = %° and 34° (other 
parameters are the same as in Fig. |S| 
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FIG. 8: Distribution of collision times of rods over the phases 
of the plate vibration for different values of the tilt angle (p. 
Solid line: approximation H{6) — (1 +cos6')/27r used for the 
theoretical analysis, where 9 is the phase of the plate velocity, 
Vpi = Vb cos 6. Parameters are the same as in Fig. |S| 



spontaneous formation of the non-uniformity of rod ar- 
rangement in an annulus at intermediate inclination an- 
gles. For small tilt, the rods are arranged in a uniform 
hexagonal-like packing (see Fig. [SJ with one rod near 
the inner wall of the gap, next near the outer wall and 
so on. This "perfect" arrangement may be perturbed by 
the presence of one or few one-rod defects due to geomet- 
rical frustration, however their cumulative effect is quite 
small and does not change the collective motion consid- 
erably. At (/) ~ 18° in most of the numerical experiments, 
the hexagonal packing spontaneously breaks and as a re- 
sult a localized region with much larger tilt emerges. At 
(j) w 25° it involves roughly half of the rods [Fig. Ele)]. 
Eventually, this region spans the whole perimeter and the 
dependence Cx (0) becomes smooth again. This defect ap- 
pears to play a significant role in slowing down the rods 
drift. This effect is exacerbated in numerical simulations 
by neglecting the rod rotation around its axis. Therefore, 
rolling of rods along the side walls is prohibited by the 
numerical model, and thus the sidewall friction is effec- 
tively enhanced. 

The second mechanism operates both in 2D and 3D, 
and it has to do with the bifurcation between long flights 
spanning more than one period of vibrations at small tilt 
angles and short flights which last one period of vibration 
at large tilt angles (see Figured). This transition occurs 
at « 30° and it leads to the noticeable difference in 
the distribution of the landing times over the phase of 
the plate vibrations for different tilt angles (see Figure 
IS)). At large angles the landing times are mostly con- 
fined to the phase interval in which the plate moves up- 
ward, whereas at smaller tilt angles there is a significant 
probability of collisions during the downward motion of 
the plate {tt/2 < 8 < 3tt/2), where 6 = mod(27r/t, 27r). 
This bifurcation explains the bump in the dependence of 
the mean vertical plate velocity at impact [see Fig. El^d)] 
and correspondingly the horizontal translation velocities 



Cx,c'^ on the tilt angle [see Fig.Efb)]. 

Overall, our simulations show that the side walls do 
play an important role in determining the magnitude of 
the horizontal velocity of rods. As seen from compar- 
ison of Figs. 0] and El the transport velocity in quasi- 
2D case is 2.5 times greater than in the annulus for the 
same values of parameters. On the other hand, they al- 
low us to develop a theoretical model of the collective 
rod motion based on the observations that, to a first ap- 
proximation, the pre- and post-collision center of mass 
translation speeds are close to each other and that the 
angular velocity before the collision can be neglected. 



IV. ROD COLLISION WITH PLATE 

In this section we derive the necessary formulas for an 
isolated collision between a rod and a vertically moving 
plate. We confine the analysis to the case of in-plane, 
eccentric, oblique frictional impact. Our derivation is 
based on the classical analysis which assumes the con- 
stant kinematic restitution coefficient (see Ref. 0). 

Consider uniform rigid rod of mass m and length I 
colliding with a plate at point O (see Fig IS)). We place 
the system of coordinates at point O so that the com- 
mon normal n coincides with the ort z and the axis of 
rod, prescribed by unit vector u, is in the x — z plane, 
u = (sin(0), 0, cos(0)). Prior to the collision, the rod has 
translational velocity (associated with the center of mass 
G) c — {cx,0,Cz), and angular velocity to = (0,a;y,0), 
and the plate has only vertical velocity V = (0,0, Vz). 
The corresponding post-coUisional velocities of rod are 
denoted by primes. 

Newton's second law gives equations for translational 
and rotational motion of the rod. In differential form 
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FIG. 9: Geometry of the collision between rod and sphere. 



they read, 



Iduj 



mdc 
I 



-u X dP 



(1) 
(2) 



where / is the moment of inertia of a rod for its rotational 
motion around the center of mass G. dP = Fdt is the 
differential of the impulse P exerted on the rod during 
the collision. 

The impulse acquired by the rod is the integral of the 
reaction force over the time of collision. The reaction 
force depends on the relative velocity at the contact point 
(CP), 



I 



V. 



(3) 



The Newton's law for the velocity of the contact point 
(CP) reads 



m— = F + a[F-u(F-u)] 



or, in projections. 



F fc2 



(4) 

(5) 
(6) 



where X, Z are coordinates of the center of mass (CM) 
of the rod, and k is the radius of gyration of the rod. 
In writing Q we assumed that the time of collision is 
so short that we can neglect the changes in the plate 
position and velocity. Depending on initial conditions, 
the impact proceeds according to one of three different 
scenarios. 

Slide. Let us denote the duration of the contact tf, 
so ^^^(0 <t <tf) > 0,^^(0) = F^{tf) = 0. We assume 
that at t — Q, Vx{Q) — uq > Q and Vz{Q) — vq < Q. 
After initial contact, the rod slides along the plate, so 
Fx = —fJ-Fz (for brevity we dropped the subscript rb of 



the friction coefficient). If 

XZ + ii{k^ + ^2) 



uo 



Pzitf) > 



(7) 



[here Pz{t) — J* Fz{t)dt] the slip in positive direction 
continues throughout the contact, and w/ = u{tf) — u^. 
While the total impulse Pz(tf) is not known a priori, it 
can be determined from the kinematic condition Vf — 
—cvq assuming that |7J) is satisfied. Then integrating 
Eq. from i = to i/ we get 



„ / N / N m,k 

P.(t/) = -Mi + ^) ,.^^.^^^^ 



and correspondingly 

Px{tf) ^Vo{l + e) 



mfik^ 



fc2 +X2 + ^iXZ' 



(8) 



(9) 



Now we can calculate the CM velocities after the contact 
using 



4 = Ca; + Px{tf)/m, 

c'z = Cz + Pzitf) /m. 



(10) 
(11) 



Replacing uq = Cx — loZ and vq = Cz — Vz + loX we get 
the center of mass (CM) velocities immediately after the 
collision 



cz - 



fi{l + e)k^(cz~Vz+u;X) 

^{k^ + Z^) + XZ 
(cz + uJx - Vz){l + e)k^ 
fi{k^ + Z^)+XZ 



(12) 
(13) 



Using (jSJ the no-stop sliding condition Q can be writ- 
ten as 

XZ + n(k^ + Z^) , , 

If the condition (|14|l is violated, at a certain time ti dur- 
ing coUision the sliding stops, Vx{ti) = 0. At t = ti, 
the horizontal CP velocity ui = and the vertical CP 
velocity is 



P +X^ + fiXZ 



(15) 



At t > ti the contact may remain at rest or reverse the 
direction of sliding. 
Slip-stick. If 



^(fc2 



Z^) > XZ, 



(16) 



the contact sticks, and the horizontal velocity u{ti < t < 
tf) = 0. During this phase, 



XZ 

k^ + Z^ 



F. 



(17) 
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Assuming again the kinematic restitution law Vf = — ewo, 
we derive the expression for the CM velocity at the end 
of the collision (see Appendix B) : 

, _ {c.^~ujZ)Z^-{c,-V,+ujX){l + e)XZ 



fc2 + X2 + Z2 



-ljZ 



(18) 

(c^ -V^+ ujX)[X^ - (fc2 + Z^)e] - (c^ - ujZ)XZ 



fc2 + X2 + Z2 



- LuX 



(19) 



Note that for the case of slip-stick, the post-collisional 
velocities are independent of the friction coefficient. 
Slip reversal. If 



fi{k^ + Z^) < xz, 



(20) 



the contact slides back after stopping. In this phase 
Fx — iiFz- Again omitting the details of derivation (see 
Appendix B) we give here the formulas for the CM ve- 
locity at the end of collision: 



Cx + {cx - ujZ) 



2^fc2(/c2 + x2 



-(cz - V,+ujX) 
c'z ^ Cz ~ {cx - ujZ) 



{fl{k^ + Z2) + XZ){k^ + X2 - /iXZ) 

/ifc2(l + e) 



k^ + X^- fiXZ 

2^iPXZ 



(21) 



-icz-Vz+u;X) 



{fl{k^ + Z2) + XZ){k^ + X2 - /iXZ) 

A:2(l + e) 



^2 - ^XZ 



(22) 



Thin rod. For a thin rod of length /, we use values 
X = lam(j)/2,Z = lcos(l)/2,k = 1/2V3, 1 = mP/l2. Let 
us first outline the boundaries of three different regimes 
(continuous slide, slip=stick, and slip reversal) in the 
((/), ji) parameter plane. While the condition for the tran- 
sition from slip-stick to slip reversal is universal, the con- 
dition of continuous slide 1)14(1 depends on the values of 
e and the ratio Cz/cx and lo. Figure 1101 shows the bi- 
furcation diagram for e = 0.8 and three different val- 
ues of —Cz/cx for uj — Q. Typically this ratio is large, 
so the regime of sliding can only be observed for small 
fi < = — Ca;[4cz(l -|- e)]~^ and either large or small (f>. 
For larger ^ > /Jc, at small the slip-stick regime occurs, 
and at larger angles there is the slip reversal regime. The 
critical angle (jjc at which the transition from slip-stick to 
slip reversal occurs, is determined from equation for (jjc 
(3 sin (/)c cos 0c = m(1 + 3 cos 02). Solving this equation, 
we obtain 




1 J9 - 16/^2 - 5f^2 
_arccos^^^ 



M^) 



(23) 



For small (j)c — f^i + 0{fi'^). Note that the critical 
angle is only dependent on the friction coefficient and 
becomes tt/2 at = 3/5. At larger fi, the slip-reversal 
scenario does not occur at any tilt angle. 



FIG. 10; Bifurcation diagrams of single rod collision with a 
plate for e = 0.8 and three different values of —Cz/cx'. 1 (a), 
2 (b), and 10 (c). 



V. COLLECTIVE MOTION OF RODS 

In order to analyze the collective motion of bouncing 
rods using the results for an individual rod collision ob- 
tained in the previous section, we have to make additional 
assumptions regarding the interactions of rods. In the 
formulation of these assumptions, we use the numerical 
and experimental observations. Referring to Figure |Bfa- 
c), we assume that in the stationary translation regime, 
Lu = 0,Cx = c'^,Cz = —c'z- Note that the latter simple 
assumption is not very accurate for large F and small 0, 
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however using it still leads to a reasonable agreement be- 
tween the theory and simulations. A more accurate set of 
closure conditions would require a detailed description of 
the complicated interactions of the rod during the flight 
between two consecutive collisions. 

Adopting these simplifications, we immediately arrive 
to the relations for the horizontal and vertical velocities 
of the rods. Note that among the three cases outlined 
above, the sliding regime cannot be realized in the regime 
of stationary translation, since it would imply a contin- 
uous decay of c^- So eventually one of the two other 
regimes will be established depending on the inclination 
angle (in a finite container, the dynamically selected in- 
clination angle is weakly dependent on the driving accel- 
eration, see inset in Fig. Il^f b)). 

Slip-stick. Assuming c'^ = Cx, c'^ = —Cz, and oj — 0, 
we get from (|18|) . (|19() in the slip-stick regime (0 < cj)c) 



, _ 2{l + e)XZV, 
" P(l-e) + 2X2 

, ^ (1 + e)PV, 

For a thin rod, Eq. (|^ yields 

, 6(1 + e) sin cos 
'^'^ ^ l-e + 6sin2 
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Slip reversal. In the slip-reversal regime i 
solve Eas. (|21|) . H22|l with constraints Cx = c'^ 
u; = 0. As a result, we get 



(24) 
(25) 

(26) 
b > (f>c), we 

Cz = -c' 



(1 + e)fiik^ + Z^)V, + (1 + e)XZV, 

(1 + e)k^V, 
k^{l-e) + 2X^ 



(27) 
(28) 



Note that the vertical velocity is again independent of 
/i and in fact coincides with (|25|l . It is easy to see that in 
the transition point from slip-stick to slip reversal regime 
where XZ — ^{k^ -\- Z'^) the values of the horizontal 
translation speed H24() and (|27|l coincide. 
For a thin rod case, we obtain from Ea. (|27|l 



(1 + e)[M(l + Scos^ (j)) + 3sin(/)C0S< 



1 



6 sin 



-14 (29) 



The vertical velocity after collision is given by 

1 + e 



1 - e + 6 sin"' 



-14 



(30) 



in both slip-stick and slip reversal regimes. Figure 1111 
shows the dependence of the normalized vertical velocity 
and the translation speed c^/14, c^/14 on the inclination 
angle 4> for the case w = 0. The transition from slip-stick 
to slip reversal dependence occurs at 4'c- 

The remaining question is what is Obviously 14 
is smaller than the amplitude Vb of the plate velocity 




40 n 60 



100 



FIG. 11: Normalized vertical velocity after collision c^/Vz 
(green dot dashed line) and the translation velocity c'^/Vz 
(black solid line) as a function of the inclination angle for 
e = 0.65, /I = 0.3. Red dashed shows the unphysical branches 
of the slip-stick and slip reversal dependencies (1181 . 12111 for 
the horizontal post-coUisional velocity. 



V = Vbcos(r2t), because the rods collide with the plate 
at different phases and not only at phase when V — Vq. 
A simple assumption which we are going to make is that 
14 = olVq with a constant fitting parameter a < 1. In 
fact our numerical simulations indicate that a is close to 
0.5 but varies slightly with (p because the landing phase 
distribution depends on (j) (see Figure|Hl), but for the sake 
of simplicity we shall ignore this dependence. The value 
a = 0.5 is obtained if we approximate the distribution 
of collision phases as H{0) = (1 -f cos6')/27r, where 9 is 
the phase of the plate velocity, Vpi = Vq cos 9. Then we 
obtain for the average plate speed at collision, 



14 = 1^0 (2^)- 



2ir 



cos 61(1 + cos 6i)d6' = yo/2 (31) 



At small < r — 1 < 1, the distribution deviates sig- 
nificantly from H{9). The rods only leave the plate for 
short flights near the top position of the plate at which 
the vertical acceleration is smaller than —g^ and the ver- 
tical velocity V is close to zero. Due to inelasticity, after 
landing the rod may perform a few more smaller bounces 
before coming to rest until the next period. While it is 
difficult to describe this regime analytically, one can ex- 
pect that 14 oc r — 1 at small F. 



VI. DISCUSSION 

In this section we compare the theoretical results with 
numerical simulations for the quasi-2D case. First, we 
tested the theoretical predictions for the isolated rod 
bounced off the plate, and found a good agreement be- 
tween formulas lO, (HHJ, (El, lEU, lES, and 
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FIG. 12: Post-coUisional CM velocities and as a function 
of tilt angle for n = 0.4 and —Cz/cx = 2. Symbols indicate 
the results of numerical simulations, and lines show theoreti- 
cal predictions for different collision scenarios using constant 
coefficient of restitution, e = 0.65. Solid lines denote slip- 
stick, dash lines denote continuous slip, and dash-dot lines 
denote slip-reversal. 



numerics (see Fig. I12() . One can clearly see the tran- 
sitions between three different regimes of rods bouncing: 
slide, slip-stick, and slip reversal. A slight difference be- 
tween the theory and simulations in the vertical veloc- 
ity at small tilt angles is related to the above-mentioned 
variations of the kinematic restitution coefficient with tilt 
angle in soft-particle MD simulations which was ignored 
in an analytical calculations. 

A comparison between the theory and MD simulations 
for the mean translation velocity is shown in Figure [T^ 
Figure ITST a') shows the horizontal translation velocity as 
a function of the mean tilt angle. The overall dependence 
is in reasonable agreement with the theory, however some 
noticeable differences are obvious. This should not be 
surprising, given the crude assumptions made to describe 
the dynamics of rods during flights between collisions. As 
mentioned above, the 'bump' visible at « 30° is related 
to a transition from the regime of 'long' flights that span 
more than one period of vibrations, to the 'short' flights 
that last a fraction of the period of vibrations. As seen in 
Figure |S| these two regimes are characterized by signif- 
icantly different distributions of the collision times over 
the vibration phases. 

Figure lT^T b) shows the F dependence of the horizontal 
translation velocity. As expected from the theory, and 
seen in experiments, the horizontal translation velocity 
is linearly proportional to F at large F. Unlike the annu- 
lar case, the translation velocity turns zero at F = 1 , and 
indeed it grows as F — 1 at small F — 1. Figure IT^ c') ad- 
dresses the question of the rod length dependence of the 
horizontal translation velocity. According to the theory 
for thin rods, Cx should be independent of L On the other 
hand, the drift should disappear when the aspect ratio 
of the rod approaches 1 (the case of spherical particles) . 
As seen in Figure ll^f c) , the translation velocity grows 



linearly at small / > 1, but this growth saturates at Z « 3 
after which the translation velocity is independent of / in 
agreement with the theory. 

There are several possible sources of discrepancies be- 
tween the theory and numerics (and experiments). First, 
in describing individual collisions we made an implicit 
assumption that the collision time is small as compared 
with the period of oscillations. This assumption may 
break for high frequency vibrations or in the regime of 
small F when rods spend a signiflcant portion of the pe- 
riod in contact with the plate. Furthermore, we used 
the simplest closure assumptions to relate the horizontal 
and vertical velocities of rods after and before the flight. 
While the relation Cx = appears to work well through- 
out the range of parameters corresponding to experimen- 
tal conditions, the other condition Cz = —c'^ holds only 
approximately. Our numerical experiments showed a sig- 
nificant (up to 50%) deviations from this simple relation 
at large F, when many inelastic collisions occur during 
the flight. We were unable to describe this effect theo- 
retically, and chose to sacrifice the accuracy of compari- 
son rather than introduce an unknown fitting parameter 

-Cz/c'z- 

Comparison between 2D and 3D simulations (Figs. 01 
I13|l shows that the characteristic translation velocity in 
3D case is 2.5 times smaller than in 2D case with the 
same material parameters. This difference may be ac- 
counted for by the frictional interaction with side walls. 
We systematically studied the dependence of the transla- 
tion velocity in the annulus on the friction coefficient Hrw 
between the rods and the side walls, and found that in- 
deed it varies strongly with firw, in particular, the trans- 
lation velocity at iirw — 0.3 is 2.5 times smaller than 
Cx at = (see Fig. I13f b)). We also analyzed the 
dependence of the translation velocity on the friction co- 
efficient with the bottom, and found that for large /i^b the 
translation velocity becomes independent of iirb- This is 
consistent with the theoretical argument that at large fj,rb 
the slip-stick scenario occurs for an arbitrary tilt angle 

As a conclusion, we studied experimentally and the- 
oretically the drift of anisotropic particles (rods) on a 
vibrated plate. The experiments in the annulus showed 
the robust drift of rods in the direction of their tilt, at the 
normalized vertical acceleration F > 1.5. For smaller val- 
ues of 1 < F < 1.5, very small reverse drift was observed. 
We developed a novel numerical algorithm which allowed 
us to study the interaction of rods in soft-particle MD 
simulations. Simulations of rods in an annulus with pa- 
rameter closely matched experiment, revealed very simi- 
lar behavior, both qualitatively and quantitatively. 

Out theoretical description of the rod translation is 
based on the detailed analysis of frictional collisions be- 
tween an individual rod and the moving plate. The ef- 
fects of collective interaction of rods during flights be- 
tween collisions are taken into account using the sim- 
plest phenomenological closure conditions based on the 
experimental findings and MD simulations. As a result. 
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closed formulas for the horizontal translation velocity are 
obtained. A direct comparison between the theory and 
experiments is complicated by the role of interaction of 
rods with frictional side walls which is unaccounted for in 
the theory. However, we found a reasonable agreement 
between the theory and numerics for quasi-2d geometry 
when rods are confined to the x — z plane with periodic 
in X boundary conditions. Since the same numerical code 
describes well the experiment in the annulus geometry, 
we infer that the theory correctly captures the mecha- 
nism of the rod translation in experiment. 



Some more subtle effects are, however, are difficult to 
model theoretically. The (very slow) reverse drift of rods 
for small F is presumably due to the small negative value 
of average Vz , however to calculate the mean Vz one needs 
a detailed knowledge of the distribution of landing times 
for small T which is difficult to obtain theoretically. 



2 3 4 5 6 
/(cm) 



FIG. 13: Results of simulations in quasi-2D geometry. A*' = 40 
rods are placed in a periodic domain of different lengths L 
which determine the tilt angle (/f>o. (a) - average horizontaf 
vefocity of rods as a function of the incfination angfe for F = 
3.3, sofid fine - theory ^ with e = 0.65, fi = 0.3; 

(6) - average horizontal vefocity of rods and average tift (inset) 
as a function F for L — 43; (c) average horizontal vefocity 
of rods and average tift as a function of the fength of rods 
I = ho + 1, where ho is the distance between the centers of 
sphericaf caps. N = 40 rods are pfaced in a periodic domain 
of length L = 41.4 [l3, the bottom is osciffating at frequency 
/ — 60Hz and acceferation is F = 2.5. 
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APPENDIX A: MOLECULAR DYNAMICS 
ALGORITHM 

In our MD algorithm, two virtual spheres of diam- 
eter d, with centers at and rj, and with velocities 
Vi and Vj, interact via normal and tangential forces, 
Fij = Fniiij+Ft, Fn = knS^^^-^nMeSvn- We introducc 
tangential spring with deflection obtained by the integra- 
tion of tangential velocity through the period of impact, 
s = J drvt, then the tangential force component is de- 
fined separately in two cases: Ft = —ktSs — jtMeSvt, for 
stick phase, and Ft — — ^rrFn^ij for slip phase. During 
the slip phase the magnitude of s is adjusted to fulfill 
|Ft| — lJ-rr\Fn\- Here = M/2 is reduced mass for 
rod-rod collision, m is the mass of the rod, 5 — d — rij 
and Vn = 'Vij • are the overlap and relative velocity 
in the direction of normal, = (r^ — rj)/ry , while tan- 
gential direction t^^ = vt/wt is specified by the relative 



of friction between rods. MD is performed in reduced 
units, all quantities are normalized by an appropriate 
combination of the diameter d, mass of virtual sphere m, 
and gravitation acceleration g. Typical values of ma- 
terial parameters are, kn = 5 10^ {mg/D),kt = 95fc„, 
= -yt = 4 10^(5/1?)^/^. The coefficients of friction for 
rod-rod and rod-bottom collisions are firr = 0.3,firb = 0.3 
respectively. Unless specified otherwise for interactions 
with walls we also used firw = 0.3. 

To expedite the integration of Newton's equation we 
used simple leapfrog algorithm [T^ with constant time 
step At = 2.0 10~^{d/gy/^ , however we tested that ap- 
plication of more accurate integration scheme such as 5*^- 
order Gear predictor-corrector did not introduce consid- 
erable changes. 

Our choice of the values of material parameters is nei- 
ther optimal for the comparison with experimental data 
nor unique. Because we observed very good agreement 
with the theoretical description for a single collision (see 
Fig. [T^ we expect our algorithm to capture details of 
short-term collision with plate. However, for a long-term 
collision it may not be accurate. 



APPENDIX B: DERIVATION OF REFLECTION 
COEFFICIENTS 

Slip-stick. The stopping condition u{ti) = gives 
the total vertical impulse exerted by a plate on a rod for 
0<t<ti, 



Pz{h)=uo 



mk 



and correspondingly 

Pxih) = -uq 



^i{k^ + Z^)+XZ 
mfik'^ 



^iik"^ + Z-^)+ XZ 
Vertical velocity at the end of the contact 
k^ + X^ + Z^ ' 



(Bl) 



(B2) 



Vf = Vi 



k^ + Z^ 



Fz = Vq 



k^+X^+ fiXZ P + X^ + Z^ 



iPz{tf)-Pz 



Assuming the kinematic restitution law Vf = —evQ and 
using l|BTjl . l|B2|l . we get 



Pzitf)^U0 



mk^ 



VQ{l + e) + 

Pxitf) = -Wo 



Hik"^ + Z"^) + X Z 
k^+X'^+ fiXZ 



^iik"^ + Z'^)+ XZ 
m^k^ 

'iKWTz^YTxz 



Uq 



tangential velocity vt 



fir 



is coefficient 



, P+X'^ + fiXZ 



(B4) 

m{k^ + Z^) 
k^+X^ + Z2 

(B5) 

mXZ 

(fc2+X2 + Z2) 
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Now we can calculate the CM velocities after the contact. 
Substituting Hlj4|l . (|lj5|l with uq = c^. — wZ and Vq = 
c,-V,+ LuX into (ITni,((ni), we get lf]l^ . (P| . 

Slip reversal. Using and the slip condition 

Fx = fiFz for <i < t < tf, we obtain the horizontal and 
vertical velocities at tf, 



Uf 



^lik^ + Z^) -XZ 



e + x^ 



(Pzitf) - P,{h)) (B6) 



Vl 



fiXZ 



{PSf)^Pz{ti)) (B7) 



Again assuming kinematic restitution condition Vf — 
—€Vq and using (|Bip . l|B2|l . we get for the impulse dur- 
ing the reversal phase 



•mkF' 
k-2 + X2 - ijXZ 

Px{tf) = - Wo 



/i(fc2 + Z2) + XZ 



(B8) 



(1 + e)vQ 



k^ + X^ + fiXZ 
fi{k'^ + Z'^)+XZ 



Uo 



^i(k'^ + Z'^) + XZ 



(B9) 



X^ - ixXZ 



(1 + e)vo + 



+ ^2 + ^LXZ 

JliWTz^Y+xz 



Uo 



The final horizontal velocity of CP 



Uf 



/i(fc2 + Z2) - XZ 

' k^+X^ - fiXZ 



k'^ + X'^ + fiXZ 

(^ + ^^"° + Mfc^ + z2) + xz ^° 

(BlO) 

Substituting ||B8l) , l|B9l) into we get the CM 

velocities after the contact for the case of slip reversal 



/ifc 



1 - 



Aifc2 



fi{k^ + Z^) + XZ 



- ujZ 



k^ +X^ - ^iXZ 
4 = Wo + Mo 



k"^ + x'^ + ^JLXZ 



Ai(fc2 + Z2)+XZ 



A;2 + X2 - ^XZ 



k^ + x^ + ^lXZ 



Substituting uq = Cx — uiZ, vq — Cz — Vz + ujX, we get 
(|aj,(|23)- 



